Nonequilibrium critical behavior of magnetic thin 
films grown in a temperature gradient 



Julian Candia and Ezequiel V. Albano ' 

1 - Instituto de Fisica de Liquidos y Sistemas Biologicos, 
(CONICET, UNLP), 59 Nro 789, 1900 La Plata, Argentina 

2 - Departamento de Fisica (UNLP), La Plata, Argentina 

July 12, 2012 

Abstract 

We investigate the irreversible growth of (2 + 1)— dimensional magnetic 
thin films under the influence of a transverse temperature gradient, which is 
maintained by thermal baths across a direction perpendicular to the direction 
of growth. Therefore, different longitudinal layers grow at different temper- 
atures between 71 and T 2 , where T x < T c hom < T 2 and T c hom = 0.69(1) is 
the critical temperature of films grown in homogeneous thermal baths. We 
find a far-from-equilibrium continuous order-disorder phase transition driven 
by the thermal bath gradient. We characterize this gradient-induced critical 
behavior by means of standard finite-size scaling procedures, which lead to 
the critical temperature T c = 0.84(2) and a new universality class consistent 
with the set of critical exponents v = 3/2, 7 = 5/2, and /3 = 1/4. In order to 
gain further insight into the effects of the temperature gradient, we also de- 
velop a bond model that captures the magnetic film's growth dynamics. Our 
findings show that the interplay of geometry and thermal bath asymmetries 
leads to growth bond flux asymmetries and the onset of transverse ordering 
effects that explain qualitatively the shift observed in the critical tempera- 
ture. The relevance of these mechanisms is further confirmed by a finite-size 
scaling analysis of the interface width, which shows that the growing sites of 
the system define a self-affine interface. 

1 Introduction 

The importance of thin film technology has been widely recognized in the realms 
of experimental and applied science, from the manufacture of electronics (layers 
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of insulators, semiconductors, and conductors from integrated circuits) to optics 
(reflective and anti-reflective coatings, self-cleaning glasses, etc) and packaging (e.g. 
aluminium-coated PET films). Indeed, the increasing role of thin films in basic and 
applied research relies on the development and refinement of nanoscale deposition 
techniques such as sputtering and molecular beam epitaxy, which allow a single layer 
of atoms to be deposited at a time [1, 2, 3, 4]. 

Since the growth temperature is one of the critical parameters in the formation 
of ordered thin films, several experiments have focused on the influence of a tem- 
perature gradient during film growth. In an early experiment by Tanaka et al. [5], 
magnetic Tb-Fe films were grown between two substrates with a temperature gra- 
dient, reporting the observation of perpendicular magnetic anisotropies and other 
gradient-driven structural features. More recently, Schwickert et al. [6] introduced 
the "temperature wedge method" where a calibrated temperature gradient of sev- 
eral hundred Kelvin was established across the substrate during co-deposition of Fe 
and Pt on MgO(OOl) and MgO(llO) substrates. These experiments generated the 
Llo ordered phase of FePt, which is currently the leading candidate material for 
ultrahigh density heat assisted magnetic recording (HAMR) and bit patterned mag- 
netic recording (BPMR) media ([7, 8] and references therein). Other experiments 
by Yongxiong et al. [9] have investigated the evolution of Fe oxide nanostructures on 
GaAs(lOO) by using a multi-technique experimental setup that included transmission 
and reflection high energy electron diffraction and scanning electron microscopy. In 
these studies, nanoscale epitaxial Fe films were grown, oxidized, and annealed using 
a gradient temperature method, which led to nanostripes with uniaxial magnetic 
anisotropy. As a result of the experimental advances on this field, many techno- 
logical applications have been envisioned as well. For instance, magneto-optical 
recording studies of signal reproduction [10] have suggested that recording media 
having multiple magnetic layers in a transverse temperature gradient may suppress 
magnetic noise from tracks adjacent to the target track during information storage 
and reproduction [11]. 

From a theoretical perspective, gradients have been studied extensively in the 
context of diffusion processes and later extended to thermal conductivity and heat 
conduction problems. The so-called gradient percolation method was originally in- 
troduced to study percolation transition models where the density is the control 
parameter [12] and later applied to a variety of problems, such as fractal diffu- 
sion fronts [13, 14, 15], overlapping disks in a concentration gradient [16], bond 
percolation for the Kagome lattice [17], invasion percolation under gravity [18], 
porous media [19], as well as in the study of vegetation distribution [20]. Very re- 
cently, the gradient method has been extended as a powerful tool to study first- and 
second-order irreversible phase transitions in far-from-equilibrium systems such as 
the Ziff-Gulari-Barshad model and forest-fire cellular automata [21, 22]. In magnetic 
systems, damage spreading processes in a temperature gradient [23] and studies of 
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several one- dimensional models [24, 25, 26, 27] have been followed by the inves- 
tigation of the kinetic Ising model in two dimensions under a variety of dynam- 
ics [28, 29, 30, 31, 32]. 

Within the broad context of these recent experimental and theoretical investiga- 
tions, the aim of this paper is to study the irreversible growth of magnetic thin films 
in a temperature gradient and to provide a full characterization of the gradient- 
induced critical phase transition. The magnetic thin film growth process under 
far-from-equilibrium conditions is investigated by using the so-called magnetic Eden 
model (MEM) [33, 34, 35], an extension of the classical Eden model [36] in which 
particles have a two-state spin as an additional degree of freedom. Earlier studies 
have shown that, growing in (d + l)-dimensional strip geometries in homogeneous 
thermal baths, MEM films are noncritical for d — 1 [34]. In contrast, for d = 2 they 
undergo an order-disorder phase transition that takes place at T^ om = 0.69(1) in the 
thermodynamic limit. The critical exponents are p hom = 1.04(16), ^ hom = 2.10(36), 
and f3 hom = 0.16(5), which intriguingly agree within error bars with the exact expo- 
nents for the kinetic Ising model [34]. Since the MEM growth process is irreversible 
and newly deposited particles are not allowed to flip and thermalize once they are 
added to the growing cluster, the observed correspondence between the MEM and 
the equilibrium Ising model remains puzzling. 

In this work, we focus on the critical case (i.e. d = 2) and show that, when 
applying a transverse temperature gradient maintained by thermal baths between 
temperatures 7\ and T 2 , where T\ < T^ om < T 2 , the system undergoes a continuous 
phase transition at a higher critical temperature (T c > T c Aom ) and with different 
critical exponents. We also develop a growth bond model and show the existence of 
bond flux asymmetries caused by the interplay of geometry and thermal bath asym- 
metries, which shed some light on the growth dynamics and explain qualitatively 
the shift observed in the critical temperature. The growth bond model analysis is 
further supported by the fact that the growing interface is self-affme, thus ensuring 
that the growing sites are correlated at all size scales. 

The rest of the paper is laid out as follows. In Section 2, we introduce the 
model and describe the Monte Carlo algorithm used to simulate MEM thin films 
in a temperature gradient. In Section 3, we present our results and a discussion. 
Finally, Section 4 consists of concluding remarks. 

2 The Model and the Monte Carlo Simulation 
Method 

The MEM in (2 + 1)— dimensions is studied in the square lattice by using a rectan- 
gular geometry L x X Ly X L z , where L z ^> L x = L y = L is the growth direction. 
The location of each spin on the lattice is specified through its coordinates (x, y, z) 
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(1 < x,y < L, 1 < z < L z ). The starting seed for the growing cluster is a plane of 
LxL parallel-oriented spins placed at z — 1 and cluster growth takes place along the 
positive longitudinal direction (i.e., z > 2). Across one of the transverse directions 
(the y— axis), a temperature gradient is applied by thermal baths at fixed tempera- 
tures linearly varying between T\ and T 2 . Therefore, in our setup each layer at fixed y 
is subjected to a constant temperature T(y) = Ti+(T2— Ti)(y— 1) / (L—l) maintained 
by a thermal bath. We adopt open boundary conditions along the y— direction, while 
continuous boundary conditions are considered along the x— direction. 

Clusters are grown by selectively adding two-state spins (S xyz = ±1) to perimeter 
sites, which are defined as the nearest-neighbor (NN) empty sites of the already 
occupied ones. Let us recall that the substrate is a 3D cubic lattice and therefore 
each lattice site in the bulk has 6 NN sites. Considering a ferromagnetic interaction 
of strength J > between NN spins, the energy E of a given configuration of spins 
is given by 

E = —— ^2 SxyzSx'y'z' , (1) 
(xyz,x'y' z') 

where the summation (xyz, x'y'z') is taken over occupied NN sites. The Boltzmann 
constant is set equal to unity throughout, and both temperature and energy are 
measured in units of J. The probability for a perimeter site at (x, y, z) to be occupied 
by a spin is proportional to the Boltzmann factor exp(— AE/T), where AE is the 
change of energy involved in the addition of the spin and T is the temperature at 
the perimeter site. 

At each step, all perimeter sites have to be considered and the probabilities of 
adding a new (either up or down) spin to each site must be evaluated. Using the 
Monte Carlo simulation method, after all probabilities are computed and normal- 
ized, the growth site and the orientation of the new spin are both simultaneously 
determined by means of a pseudo-random number. Notice that the MEM's growth 
rules require updating the probabilities at each time step and lead to very slow al- 
gorithms compared with analogous equilibrium spin models. Since the observables 
of interest (e.g. the mean transverse magnetization along the x— direction and its 
higher moments) require the growth of samples with a large number of transverse 
planes of size LxL, clusters having up to 10 9 spins have typically been grown for 
lattice sizes in the range 12 < L < 96. Also, let us point out again that, although 
Eq. (1) resembles the Ising Hamiltonian, the MEM is a nonequilibrium model in 
which new spins are continuously added, while older spins remain frozen and are 
not allowed to flip, detach, nor diffuse. 

As in the case of the classical Eden model, the magnetic Eden model leads to a 
compact bulk and a self-affine growth interface [33] (see Sect. 3.4 for a detailed finite- 
size scaling analysis of the interface width). The growth front may temporarily create 
voids within the bulk, usually not far from the rough growth interface. However, 
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Figure V. Snapshot showing a longitudinal slice for a fixed value of the transverse 
coordinate x. A temperature gradient between T\ = 0.5 (bottom) and T 2 = 1.5 (top) 
is maintained along the transverse axis y. The system grows along the longitudinal 
z > direction in a semi-infinite (2 + l)-dimensional strip substrate. Red (black) 
sites represent up (down) spins, while empty sites are shown in white. 

since the boundaries of these voids are also perimeter sites, they ultimately become 
filled at some point during the growth process. Hence, far behind the active growth 
interface, the system is compact and frozen, and the different quantities of interest 
can thus be measured on defect-free transverse planes. 

Notice that the growth of magnetic Eden aggregates in (2 + l)-strip geometries 
is characterized by an initial transient length ~ L 2 (measured along the growth 
direction, i.e. the z— axis) followed by a nonequilibrium stationary state that is 
independent of the initial configuration [34] . We considered starting seeds formed by 
L x L up spins (i.e. S xy \ = 1) but any choice for the seed leads to the same stationary 
states for z 3> It- By disregarding the transient region, all results reported in this 
paper are obtained under stationary conditions. Also notice that, since the films are 
effectively semi-infinite and the substrate length along the growth direction plays 
no role, the only characteristic length is the transverse linear size L. 

3 Results and Discussion 

3.1 Gradient-Driven Continuous Pseudo-Phase Transitions 
in Finite-Size Films 

Let us begin by considering a fixed gradient between temperatures T\ = 0.5 and 
T2 = 1.5. The effect of changing this gradient will be discussed later. Figure 1 shows 
the snapshot of a longitudinal slice for a fixed value of the transverse coordinate x. 
The temperature gradient is applied along the transverse y— direction, while the 
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Figure 2: Symmetrized magnetization probability distributions for a system of linear 
size L = 64 and different layer temperatures between Ti = 0.5 and T 2 = 1.5. The 
sharp peaks near m ~ ±1 for T = 0.80 have been truncated. 

system grows from left to right along the longitudinal z > direction. Notice that 
the bulk grows compact, because although voids and holes in the bulk may eventually 
occur, they ultimately become filled at some point during the growth process. The 
bottom layers grow in contact with thermal baths at cold temperatures, which favor 
the formation of well-ordered spin domains. In contrast, the top layers grow in 
contact with hot thermal baths that promote bulk disorder. As will be shown below, 
the interplay of model growth dynamics, geometry, and thermal bath asymmetries 
lead to the onset of gradient- driven order-disorder critical phase transitions that can 
be quantitatively characterized. 

In order to take into account the asymmetries introduced by the temperature 
gradient, we can quantify the degree of order in the system by considering the 
magnetization of transverse columns at constant temperature (i.e. along the a;- axis): 

m (y) = tE S *yz ■ ( 2 ) 

Figure 2 shows the probability distributions of m for a system of linear size 
L = 64 growing in a temperature gradient between 7\ = 0.5 and T 2 = 1.5. Notice 
that different plots correspond to different layers and, therefore, to different tem- 
peratures within the gradient's range, as indicated. As expected for a continuous 
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Figure 3: Mean absolute magnetization as a function of the layer temperature for 
different system sizes, as indicated. 

order-disorder phase transition, the low-temperature distributions are bimodal and 
peaked at the spontaneous magnetization m = ±m sp (0 < m sp < 1). As the temper- 
ature is increased, the peaks approach each other and merge smoothly, ultimately 
leading to a Gaussian distribution peaked at m = for high temperatures, which is 
characteristic of the disordered phase. Indeed, the smooth shift of the distribution 
maxima across T ~ T C (L), from the low-temperature nonzero spontaneous magne- 
tization m = ±m sp to the high-temperature Gaussian centered at m = 0, is the 
signature of true thermally-driven continuous phase transitions [37]. 

Notice that the distributions in Figure 2 are symmetrical, since there exists a 
finite probability for fluctuations to grow and switch the magnetization from m ~ 
+m sp tom~ — m sp and viceversa. Since Monte Carlo simulations are restricted to 
finite samples, the standard procedure to avoid these shortcomings due to finite-size 
effects is to average the absolute value of the order parameter [38]. In this context, 
the appropriate order parameter is the mean absolute magnetization of transverse 
columns at constant temperature, i.e. 

(HKy) = (^\T / S xyz \) z , (3) 

where (...) 2 denotes averages along the growth direction z > within the stationary 
region. Figure 3 shows plots of the mean absolute magnetization as a function of 
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Figure 4: Susceptibility as a function of the layer temperature for different system 
sizes in the range 12 < L < 96. As expected for a critical system, the peaks become 
sharper and higher as L is increased. 

the layer temperature for different system sizes in the range 12 < L < 96. For 
any given system size, at low temperatures the system grows ordered and the mag- 
netization is close to unity, while at higher temperatures the disorder sets on and 
the magnetization becomes significantly smaller. However, fluctuations due to the 
finite system size prevent the magnetization from becoming strictly zero above the 
critical temperature, so the transition between the low-temperature ordered phase 
and the high-temperature disordered phase becomes smoothed out and rounded. As 
expected, larger systems are less affected by finite-size effects and display sharper 
transitions. 

Strictly speaking, however, these results just show evidence of pseudo-phase tran- 
sitions, which might be precursors of true phase transitions taking place in the 
thermodynamic limit. In the following, we will characterize in more detail this 
pseudo-critical phenomenon by measuring other observables on finite-size systems. 
In the next Subsection, we will use standard finite-size scaling procedures to estab- 
lish the existence of a non-trivial critical temperature in the L — > oo limit, as well 
as to calculate critical exponents that describe the behavior of the infinite system 
at criticality. 
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Let us now consider the magnetic susceptibility x, given by 

x = y (> 2 > - <H> 2 ) • (4) 

For equilibrium systems, the susceptibility is related to order parameter fluctua- 
tions by the fluctuation-dissipation theorem. Although the validity of a fluctuation- 
dissipation relation in the case of a nonequilibrium system is less evident, we will 
assume Eq. (4) to hold also for the MEM. Indeed, as shown in earlier studies of 
nonequilibrium spin models [39, 40], this definition of x proves very useful for ex- 
ploring critical phenomena under far-from-equilibrium conditions. In Section 3.2, we 
will characterize the critical behavior in the thermodynamic limit through critical 
exponents and finite-size scaling relations by applying the equilibrium theory to our 
far-from-equilibrium model. 

Figure 4 shows plots of x vs T for different system sizes, as indicated. As with 
the thermal dependence of the order parameter shown in Figure 3, the peaks of the 
susceptibility become rounded and shifted, indicating the existence of pseudo-phase 
transitions in finite-size MEM thin films. Increasing the system size, the peaks 
become sharper and higher, as expected for a critical system. 

Since the results presented thus far considered a fixed gradient between the 
temperatures 7\ = 0.5 and T 2 = 1.5, let us now investigate the effects of changing 
the gradient span AT = T 2 — T\. With this aim, we keep the same temperature 
for the thermal bath at the cold end (7\ = 0.5) and consider a substantially higher 
temperature for the thermal bath at the hot end (T 2 = 2.5). 

Figure 5 compares the mean absolute magnetization for these two different gra- 
dient ranges (i.e. AT = 1,2) for systems of size L — 12 and L = 96. We observe 
that increasing the gradient span shifts the magnetization profiles towards higher 
temperatures. That is, the temperature of a given layer does not uniquely deter- 
mine its degree of order, since the mean magnetization of the layer also depends on 
the overall gradient span under which the film grows. Similar shifts towards higher 
temperatures are also seen in higher-order moments of the order parameter prob- 
ability distributions, such as the susceptibility and Binder's fourth-order cumulant 
(not shown here for the sake of space). Alternatively, we can compare systems of 
different sizes and gradient spans such that the local gradients 5 = AT/L are the 
same. The inset to Figure 5 shows a comparison between a system of size L — 16 
and gradient span AT = 1 (dashed lines) and another system of size L = 32 and 
gradient span AT = 2 (dotted lines), both of which have the same local gradient 
5 = 1/16. The solid line corresponds to the mixed case L = 32 and AT = 1 (i.e. 
5 = 1/32). We observe that the systems with the same local gradient have simi- 
lar magnetization in layers at intermediate temperatures (i.e. approximately in the 
range 0.8 < T < 1.3). However, the magnetization profiles for equal-5 systems differ 
noticeably in layers closer to the borders of the sample. 
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Figure 5: Comparison of magnetization profiles vs T for different gradient 

spans (AT = 1,2) and system sizes (L = 12,96). The arrows indicate the corre- 
sponding shifts in the finite-size critical temperature T C (L). Inset: Magnetization 
profiles for systems with the same local gradient 5 = 1/16, namely L = 16, AT = 1 
(dashed line) and L = 32, AT = 2 (dotted line). For comparison, the mixed case 
L = 32, AT = 1 (solid line) is also shown. 

The arrows in Figure 5 show that the shifts for smaller systems are signifi- 
cantly larger than the corresponding shifts for larger systems. Defining the finite-size 
pseudo-critical temperature T C (L) as the temperature corresponding to (|m|) = 0.5, 
the shift for L = 12 is AT C = 0.18, while the corresponding shift for L = 96 is 
AT C = 0.07. This observation suggests that differences arising from changing the 
gradient span might just reflect finite-size effects that vanish in the L — > oo limit. 
In the next Subsection, we study the critical behavior of MEM films and confirm 
that, in fact, these differences are merely finite-size effects that become irrelevant in 
the thermodynamic limit. 

3.2 Characterization of Gradient-Driven Critical Behavior 
in the Thermodynamic Limit 

So far, we have found evidence for the existence of a gradient- driven order-disorder 
phase transition from the analysis of order parameter probability distributions, the 
order parameter mean absolute value (magnetization) and its fluctuations (suscep- 
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Figure 6: Binder cumulant as a function of the layer temperature for different system 
sizes, as indicated. The inset shows the cumulant intersections for the larger systems 
(32 < L < 96), which determine T c = 0.84(2). 

tibility) in the growth of finite-size magnetic films. In this Subsection, we apply 
standard finite-size scaling techniques to show the existence of this phase transition 
in the thermodynamic limit and to determine critical exponents that characterize 
the system's critical behavior and universality class. 
The Binder cumulant, defined by 

6{m z ) z 

is a fourth-order cumulant dependent on the variance and the kurtosis of the order 
parameter probability distribution. One important property of the Binder cumu- 
lant is that, for large system sizes, the low-temperature, ordered region tends to 
the value 2/3, while the high-temperature, disordered region tends to 0. Thus, in 
the thermodynamic limit, the function becomes discontinuous exactly at the critical 
temperature. Moreover, since for second-order phase transitions, the scaling prefac- 
tor of the cumulant is independent of the sample size, plots of U4 versus the control 
parameter lead to a common (size-independent) intersection point that corresponds 
to the location of the critical value of the order parameter in the thermodynamic 
limit [41]. 

Figure 6 shows the Binder cumulant as a function of the layer temperature for 
different system sizes in the range 12 < L < 96. The Inset to Figure 6 shows 
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a detailed view of the same data for the largest lattices (32 < L < 96), where the 
intersection region is indicated by a grey vertical strip. Based on this observation, we 
determine the critical temperature in the L — > oo limit as T c = 0.84(2). Interestingly, 
this value is significantly higher than the corresponding critical temperature for the 
MEM growing in an homogeneous thermal bath (i.e. in the absence of a temperature 
gradient), namely T^ om = 0.69(1) [34]. In the next Subsection, we will explore the 
growth dynamics and explain qualitatively this shift in the critical temperature as 
due to ordering effects caused by a net transverse growth bond flux induced by 
thermal asymmetries. 

Notice also that, by fixing the temperature range (i.e. the temperatures T\ 
and T-i) and increasing L, we are effectively considering different gradients 5 = 
(T 2 — Ti)/ L that become smaller as L is increased. Figure 6, which shows a fixed point 
in the Binder cumulants as the gradients are changed, provides therefore quantitative 
evidence for the existence of a gradient-independent phase transition taking place 
at the temperature T c . 

According to the finite-size scaling theory, developed for the treatment of finite- 
size effects at criticality under equilibrium conditions [42, 43], the difference between 
the true critical temperature, T c , and the effective pseudo-critical one, T C (L), is given 
by 

\T c -T c (L)\<xL- 1 ' 1 ', (6) 

where v is the exponent that characterizes the divergence of the correlation length at 
criticality. As mentioned above, we define the finite-size pseudo-critical temperature 
T C (L) as the temperature corresponding to (\m\) = 0.5. 

Let us point out that, given the lack of a comprehensive theory of non-equilibrium 
phase transitions, concepts and definitions developed in the context of equilibrium 
phenomena are customarily borrowed and applied to far-from-equilibrium phenom- 
ena as well. For a review of standard methods, see e.g. [44, 45, 46]. Indeed, although 
this approach is ad-hoc and lacks the theoretical foundations of equilibrium systems, 
it has been used extensively in the literature and has become a powerful means of 
advancing our knowledge within the realm of non-equilibrium phenomena. For in- 
stance, numerical methods such as Monte Carlo simulations or series expansions 
are restricted to finite systems and it is therefore important to understand how 
far finite-size effects influence the properties of the system. As known from equi- 
librium statistical mechanics, finite-size effects are particularly strong close to the 
critical point, where the spatial correlation length becomes comparable with the 
linear dimensions of the system. By introducing the system size as an additional 
parameter, finite-size scaling laws are used to characterize the steady state of finite 
far-from-equilibrium systems through appropriate scaling exponents, such as, for in- 
stance, the exponent v in Eq.(6) above. Moreover, this procedure allows us to define 
universality classes of non-equilibrium systems, as reviewed e.g. in Refs. [47, 48, 49]. 

Figure 7 shows log-log plots of the finite-size pseudo-critical temperatures T C (L) 
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Figure 7: Log- log plots of the finite-size pseudo-critical temperatures T C (L) as a 
function of the inverse of the system linear size, T _1 , for AT = 1 (filled circles) 
and AT = 2 (open circles). By means of standard finite-size scaling analysis, we 
obtain the critical exponent v = 1.53(6) (see more details in the text). Inset: plot of 
(|m|) x I?l v vs |T - T c | x L l / V (with f3 = 0.26) showing a data collapse for AT = 1 
and different system sizes in the range 32 < L < 96. 

as a function of the inverse of the system linear size, L~ l , for different gradients and 
system sizes, as indicated. By rewriting the finite-size scaling relation as 



we performed different least-squares fits to the data using the mean, upper-bound 
and lower-bound values for the critical temperature in the thermodynamic limit. 
The nonlinear least-squares fitting procedure was implemented using the Levenberg- 
Marquardt minimization method [50] and the results from each independent fit are 
reported in Table 1. The errors in the table are determined by the fitting algorithm 
and take into account the statistical errors for each datapoint. Figure 7 shows that 
the finite-size scaling relation fits the data very well within error bars for the range of 
values for the critical temperature that was derived from the intersection of Binder's 
cumulants. From these fits, we obtain the critical exponent v = 1.53(6), where the 
error bars reported reflect the errors derived from the evaluation of T c as well as the 
statistical errors. Notice that the data for different gradients tends to converge in 
the L — > oo limit, therefore confirming that differences arising from changing the 



T C (L) =T C + Ax L- 1 '", 



(7) 
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0.82 


2.86(3) 


1.54(2) 


1 


0.84 


2.90(3) 


1.49(2) 


1 


0.86 


2.78(3) 


1.48(1) 


2 


0.82 


3.70(3) 


1.58(1) 


2 


0.84 


3.65(3) 


1.56(2) 


2 


0.86 


3.66(4) 


1.53(2) 



Table 1: Results from fitting the data to the finite-size scaling relation, Eq.(7). 
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Figure 8: Log- log plots of the susceptibility maxima as a function of the system 
linear size for AT = 1 (filled circles) and AT = 2 (open squares), where statistical 
errors for each datapoint are smaller than the symbol size. The solid lines are finite- 
size scaling fits that yield j/u — 1.66(3). The insets display plots of xL~"'^ u vs 
|T - T c \L l ' u (for the AT = 1 case) showing separately the data collapse for (a) the 
low-temperature branch and (b) the high-temperature branch. 



gradient are finite-size effects. On the other hand, finite-size scaling theory predicts 
that plots of i\m\)LP/ v vs \T-T c \L l l v for different lattice sizes should collapse near 
the critical region. The inset to Figure 7 shows the data collapse obtained by using 
(3 = 0.26 (that is determined from the hyperscaling relation, as explained below) 
with two separate branches corresponding to the low- and high-temperature regions. 
An additional characterization of the critical behavior of this system can be 
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obtained by calculating the critical exponent 7, which describes the divergence of the 
susceptibility at the critical point. Using again the finite-size scaling theory [42, 43], 
the exponent ratio 7/1/ is related to the peak of the susceptibility measured in finite 
samples of size L by 

W oc W' v . (8) 

The symbols in Figure 8 correspond to the maxima of x plotted against the 
system linear size for different gradients, as indicated, while the solid lines are fits to 
the data using the scaling relation from Eq.(8). Statistical errors for each datapoint 
are smaller than the symbol size in the Figure. The fitting procedure (which, as 
mentioned above, was implemented as a nonlinear least-squares algorithm using the 
Levenberg-Marquardt minimization method [50]) yields 7/1/ = 1.66(3), where the 
error bars reflect the statistical errors from the fits. Using this ratio and the value 
already obtained for u, we determine 7 = 2.54(11). The insets to Figure 8 display 
plots of xL~ l/v vs \T - T c \L l / u for AT = 1 and different lattice sizes in the range 
32 < L < 96. Using the critical temperature as determined by the susceptibility 
peaks, the data collapse is shown separately for (a) the low-temperature branch and 
(b) the high-temperature branch. In the former case, data from low-temperature 
layers near T\ = 0.5 depart from the collapse and have been removed. However, 
the collapse near the critical region is remarkable and agrees very well with the 
expectations from the finite-size scaling theory. By replacing the exponents v and 7 
in the hyperscaling relation dv — 2(3 — 7 = with d = 2, we determine the exponent 
(3 = 0.26(8), where the error is determined from standard error propagation applied 
to the hyperscaling relation. Recall that we anticipated this value of (3 when we 
considered the data collapse of the scaled magnetization (see the inset to Figure 7 
above). The excellent data collapse near the critical region confirms the consistency 
and robustness of the obtained results. 

As a summary, Binder's cumulant method and finite-size scaling analysis allowed 
us to characterize quantitatively the critical behavior of nonequilibrium magnetic 
films growing in a temperature gradient. We found that differences arising from 
changing the gradient are due to finite-size effects that vanish in the thermody- 
namic limit. The system's critical temperature is T c = 0.84(2), significantly higher 
than the critical temperature for films grown in an homogeneous thermal bath, 
T^ om = 0.69(1) [34]. The critical exponents are v = 1.53(6), 7 = 2.54(11), and 
[3 = 0.26(8). Based on our findings, we conjecture that magnetic Eden films grow- 
ing in a temperature gradient belong to a new universality class characterized by 
critical exponents v — 3/2, 7 = 5/2, and (3 = 1/4. In contrast, the critical exponents 
for magnetic Eden films grown in an homogeneous bath agree within error bars with 
the exact exponents for the Ising model in d = 2 [34], namely v — 1,7 = 7/4, and 
= 1/8. 
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Figure 9: Mean growth bond flux components as a function of the gradient span 
AT for L = 32 and T x = 0.5. Asymmetries due to the temperature gradient and the 
substrate geometry are responsible for net bond fluxes along the y and z directions. 
Along the transverse direction x the system is fully symmetric, so no net bond fluxes 
are observed regardless of AT. 

3.3 Growth Bond Model and Bond Flux Asymmetries 

In this Subsection, we explore the growth dynamics by means of a simple bond 
representation. Let us recall that the MEM's growth process adds new spins, which 
are deposited one by one to the growing cluster. Although voids and holes may form 
within the bulk, ultimately all sites become filled. Hence, to each pair of neighboring 
sites, we can assign a directed bond that points from the earlier occupied site to the 
later occupied one. The components of the bond flux field at a site (x, y, z) are 



where b[s\, s?\ = 1 if the bond points from Si to S2, and b[s\, S2] = —1 if the bond 
points from S2 to s±. 

Figure 9 shows the x—, y—, and z— components of the mean bond flux (0) as a 
function of the gradient span AT for L = 32 and 7\ = 0.5. As expected from the 
symmetry along the transverse x— direction, there is no net bond flux in x: (<p x ) = 



defined as: 



<f>x(x,y,z) 
(j) y (x,y,z) 
0z(x,y,z) 



b[(x,y,z), (x + l,y,z)] , 
b[(x,y,z), (x,y+ l,z)] , 
b[(x,y,z),(x,y,z + l)\ , 



(9) 
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regardless of AT. For AT = 0, the system is also symmetric along y, so no net bond 
flux is observed. When a gradient is applied, however, this symmetry is broken. 
Since the growth probabilities depend on the Boltzmann factor exp(— AE/T(y)), 
where T(y) is the layer's temperature, the thermal asymmetries introduced by the 
gradient favor spin deposition on the colder layers. This phenomenon is captured 
by the observed net bond flux ((fr y ) > 0. Indeed, as shown in Figure 9, the thermal 
asymmetries cause ((fry) to grow steeply up to ((fr y ) ma x ~ 0.75 followed by a moder- 
ate decrease for larger gradients, which is due to the onset of bulk disorder within 
the hotter layers. Since the net transverse growth bond flux is directed from the 
ordered (cold) layers towards the disordered (hot) ones, this gradient-induced trans- 
verse ordering mechanism causes the system's critical temperature to increase from 
T c hom = 0.69(1) to T c = 0.84(2). On the other hand, for AT = 0, (<fr z ) = 1 due to 
the longitudinal asymmetries in the substrate (i.e., the semi-infinite strip geometry 
constrains the system to grow along the z > direction). However, when the trans- 
verse gradient is applied, two effects contribute to decrease ((fr z ): (i) the onset of the 
transverse bond flux, which creates transverse domains in the active perimeter and 
causes some of the added spins to grow backwards; (ii) the bulk disorder induced 
in the hotter layers (which also causes ((fry) to decrease, as discussed above). 

The mean fluxes shown in Figure 9 were averaged over sites at different temper- 
atures. In order to gain further insight, let us now investigate the dependence of the 
bond fluxes on the layer temperature T(y) = Ti + (y — 1) x (T2 — Ti)/(L — 1) (where 
1 < y < L). Since we are considering different gradient spans for a fixed system size, 
it is actually more convenient in this case to plot the bond fluxes as a function of the 
transverse coordinate y. We already observed that the transverse bond flux along x 
is null regardless of temperature, so we will focus on the growth bond fluxes along 
the y— and 2— directions. Figure 10(a) shows the upwards bond flux ((fr y (y)) vs y 
for L = 32, Ti = 0.5, and different gradient spans AT, as indicated. In the absence 
of a gradient, the flux is directed downwards at the bottom and upwards at the 
top, yielding zero net bond flux. Indeed, because of the open boundary conditions, 
empty perimeter sites at the confinement walls experience a missing-neighbor effect 
and the system grows preferentially along the center of the film as compared to the 
walls. When a gradient is applied, the flux grows steeply in the upwards direction, 
as expected. However, for larger gradients, the hotter thermal baths are capable 
of inducing disorder in the bulk and partially break the upwards bond flux on the 
upper layers. Indeed, this phenomenon causes the overall bond flux (<p y ) (averaged 
over all layers) to decrease for large values of AT, as discussed above. Similarly, 
Figure 10(b) shows the forward bond flux ((fr z (y)) vs y for L = 32, T\ = 0.5, and 
different values of AT, as indicated. For AT = 0, there is a slight missing-neighbor 
effect for the sites near the confinement walls. Since forward growth is mostly driven 
by the substrate asymmetry, this effect is much less noticeable that in the flux along 
y. When the gradient is applied, two effects contribute to reduce the longitudinal 
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Figure 10: Bond flux components as a function of the transverse coordinate y/L 
for L = 32, Ti = 0.5, and different gradient spans, as indicated: (a) transverse flux 
along y; (b) longitudinal flux along z. 

flux, as discussed above. One of them, which is dominant for small gradients, is 
due to the formation of transverse domains along y, causing some backwards depo- 
sition when the bulk is filled in. The other mechanism, which is dominant for larger 
gradients, is due to the onset of bulk disorder in the hotter layers. 

Previously, we pointed out the fact that the gradient- driven order-disorder phase 
transition occurs at a temperature that is significantly higher than the corresponding 
critical temperature for the MEM growing in an homogeneous thermal bath (i.e. in 
the absence of a temperature gradient). The results presented in this Subsection 
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allow a qualitative explanation for the shift in the critical temperature. Indeed, 
the temperature gradient breaks the transverse symmetry along the y— direction, 
causing the onset of a net transverse growth bond flux. This flux is directed from 
the ordered layers (that grow at T < T c ) towards the disordered layers (that grow at 
T > T c ), thus expanding the low-temperature ordering effects across layers at higher 
temperatures. This gradient- induced transverse ordering mechanism increases the 
system's critical temperature. 

3.4 Scaling behavior of the growth interface 

According to the analysis presented in Section 3.3, the system's critical temperature 
is increased due to gradient-induced transverse ordering mechanisms that originate 
in the cold layers. It could be argued, however, that, since the distance between 
the cold layer at 7\ and a layer at a fixed temperature T > 7\ becomes larger as 
L is increased (while keeping fixed the gradient span AT), then the transverse flux 
effects may become weaker and eventually negligible in the large- L limit. In this 
Subsection we show that the growth interface is self-affme and its shape is stable 
and independent of size. Therefore, we confirm that the effects described in Section 
3.3 are still relevant in the thermodynamic limit. 

In order to track the evolution of the growth interface, we compute the position 
of the perimeter sites in the active region every time a monolayer of L 2 spins is 
deposited. The interface width at time t is defined by 



where the sum is taken over the N p perimeter sites in the active growth region, Zi is 
the longitudinal coordinate of the i—th perimeter site, and z c = (1/N p ) J2i z i is the 
center of the interface. Figure 11 shows the interface width as a function of time 
for different lattice sizes in the range 12 < L < 96, where the time unit corresponds 
to the deposition of L 2 spins. After a short transient period, the interface reaches a 
stable saturation width, w sat , analogously to other surface growth phenomena [51, 
52]. The Inset to Figure 11 shows the dependence of w sat on the system size L, where 
the statistical errors are smaller than the size of the symbols. The fit to the data 
(solid line) shows that w sa t oc L a , where the roughness exponent is a = 1.01 ± 0.01. 
That is, the saturation width scales linearly with the system size. 

By subtracting the interface center, z c , we can compute the average interface 
profile in the stationary regime, as shown in Figure 12. We find that, when scaled 
by the mean interface width, interface profiles for different system sizes collapse 
into a universal shape for the growth interface. This nearly-linear universal shape is 
qualitatively consistent with the roughness exponent a — 1, as it has quantitatively 




(10) 
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Figure 11: Interface width as a function of time, where the time unit is the deposition 
time of L 2 spins. Inset: log-log plot of the saturation width, w sat , as a function of 
system size. The statistical errors are smaller than the size of the symbols. The 
fit to the data is shown by the solid line and agrees with a linear scaling w sat oc L 
(more details in the text). 

been determined. In passing, notice also that this universal profile shows a good 
agreement with the instantaneous snapshot displayed in Figure 1. Thus, we con- 
clude that the active growth interface is self-affme and has universal features: the 
detailed analysis of bond flux asymmetries presented in Section 3.3 for a fixed lattice 
size (L = 32) remains valid for larger systems. In particular, our analysis focused 
on the influence exerted by the low-temperature layers into higher-temperature lay- 
ers, which therefore is a gradient-induced growth mechanism still relevant in the 
thermodynamic limit. 

In addition to the roughness exponent a, self-afnne interfaces are also charac- 
terized by a fractal dimension df. In fact, within short length-scales such that 
Az <C £, where Az is the longitudinal interface distance along the growth direction 
between two points separated by a transverse distance £, the fractal dimension is 
df — 2 — a [52]. In the long length-scale limit, moreover, the fractal dimension of 
the self-affme interface is df = 1, irrespective of its roughness [52]. Therefore, we 
conclude that the growth interface of the MEM in a thermal gradient is df — 1 at 
both short and long length-scales. Indeed, these conclusions are in full agreement 
with the nearly-linear shape of the growing interface that is consistent with a unitary 
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Figure 12: Collapse of the scaled average growth profile in the stationary regime for 
different lattice sizes, as indicated. 

fractal dimension at all length-scales, as seen in Figure 12. 

Here, it is useful to compare the obtained results with those from related, well- 
studied growth models. For the standard MEM growing in a homogeneous tempera- 
ture bath at high temperatures, the attachment of spins is a stochastic (random) pro- 
cess that becomes independent of the interaction energy and the temperature [33]. 
Thus, in this limit, the growth interface of the standard MEM [33, 34, 35] agrees 
with that of the classic Eden model [36], which is well known to belong to the 
Kardar-Parisi-Zhang (KPZ) universality class [52]. The most accurate simulation 
results for the KPZ model in (2 + 1)— dimensions yield a = 0.393 ±0.003 [53], which 
agrees well with some of the formerly reported values for KPZ [54] and the Eden 
model [55]. Using this value for the roughness, the fractal dimension of the growth 
interface of the MEM in an homogeneous thermal bath at high temperatures crosses 
over from df ~ 1.6 (at short length-scales) to df — 1 (at long length-scales). More- 
over, we can safely expect that the self-affine properties of the growth interface of 
the standard MEM (in an homogeneous bath) be independent of the temperature 
within a wide range around and below the critical temperature, at least insofar the 
occurrence of a layering/roughening transition, in the sense of that observed in the 
3— dimensional Ising model [56, 57], can be neglected. Therefore, we conclude that 
the fractal and self-affine characteristics of the growth interface of the MEM in a 
constant temperature bath are quite different than those of the same model growing 
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in a temperature gradient. 

On the other hand, due to the fact that the MEM grown under a temperature- 
gradient constraint exhibits a second-order transition, one may also consider the 
self-affine properties of the interface between the ordered and the disordered phases. 
Although for systems under equilibrium conditions a useful (alternative) approach 
is the evaluation of the damaged interface [23], the damage spreading technique can 
not straightforwardly be applied for the evaluation of an interface in an irreversible 
growth model. Furthermore, the implementation and application of a cluster count- 
ing algorithm [12, 15, 16, 18, 19, 21, 22, 58, 59] to our model would require formidable 
computational task that is beyond the aim of this paper. Additional shortcomings 
for this kind of calculations are the definition of the suitable cluster (e.g. Swendsen- 
Wang vs physical clusters) and the occurrence of noticeable corrections to scaling 
that one needs to evaluate in order to obtain reliable exponents [15], which is also a 
task that lies beyond our computational capabilities. However, from heuristic argu- 
ments based on the standard scaling relationship a or d~dis — + v) [12, 15, 23], 
where a or d-dis is the roughness exponent of the order-disorder interface, we can con- 
jecture that a or d-dis = 3/5, which yields a self-affine order-disorder interface with 
short length-scale fractal dimension dj d ~ dis = 7/5. For comparison, by applying the 
same scaling relationship to the standard MEM growing in an homogeneous thermal 
bath, we obtain a ord -dis = 0.51 ±0.09 and dj d ~ dis = 1.49±0.09. Thus, although the 
growth interface is very significantly affected by the temperature gradient compared 
to the thermally homogeneous system, the geometry of the order-disorder interface is 
not so markedly affected by the gradient constraint and the results for both systems 
agree within error bars. 

4 Conclusions 

In this work, we studied magnetic thin films growing under far-from-equilibrium 
conditions in (2 + l)-dimensional strip geometries, where a temperature gradient is 
applied across one of the transverse directions. We modeled the thin film growth 
process by means of extensive Monte Carlo simulations performed on the magnetic 
Eden model (MEM), in which spins are deposited on a growing cluster with proba- 
bilities dependent on a ferromagnetic, Ising-like configuration energy. 

Firstly, we studied the thermal dependence of order parameter probability dis- 
tributions, the order parameter mean absolute value (magnetization), the order 
parameter fluctuations (susceptibility) and its higher moments (Binder cumulant) 
on finite-size magnetic films, which showed the existence of gradient-driven pseudo- 
phase transitions. Secondly, we applied Binder's cumulant method and finite-size 
scaling analysis in order to characterize quantitatively the critical behavior of MEM 
films growing in a temperature gradient. The system's critical temperature is 
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T c = 0.84(2), significantly higher than the MEM's critical temperature when grow- 
ing in an homogeneous thermal bath, namely Tj} om = 0.69(1) [34]. The critical 
exponents are v = 1.53(6), 7 = 2.54(11), and f3 = 0.26(8), which also differ from 
the MEM's exponents in the absence of a temperature gradient [34]. By changing 
the gradient span, we observed finite-size effects that vanish in the thermodynamic 
limit. Hence, the critical temperature and exponents are universal for MEM films 
growing in a temperature gradient. We also investigated the system's growth dy- 
namics by means of a bond model. We found that the interplay of geometry and 
thermal bath asymmetries leads to growth bond flux asymmetries and the onset of 
transverse ordering effects that explain qualitatively the shift observed in the critical 
temperature. Finally, we analyzed the self-affine growth interface and obtained a 
collapse of the scaled average growth profile in the stationary regime for different 
lattice sizes, which shows that growth bond flux asymmetries play a relevant role in 
the model's growth dynamics even in the thermodynamic limit. 

In the context of a great experimental and theoretical interest in magnetic sys- 
tems growing in temperature gradients, as well as a wide variety of technological 
applications that benefit from these efforts, we hope that this work will contribute 
to the progress of this research field and stimulate further work. 
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